#ifndef CUBIC_PP_SPINES_H
#define CUBIC_PP_SPINES_H
#include <iostream>

#include "Spline.h"
//完全PP样条
class CubicPPSplinescomplete:public PPSplines{
    private:
    double m_1;//左端点一阶导
    double m_n;//右端点一阶导
    std::vector<double> m;
    public:
    CubicPPSplinescomplete(const std::vector<double>& x_0,const std::vector<double>& f_0,const double& m_1,const double& m_n):PPSplines(x_0,f_0),m_1(m_1),m_n(m_n){
        CalculateCoefficients();
    }
    void CalculateCoefficients()override;
    std::vector<double> getFirstDerivatives(){ 
        CalculateCoefficients();
        return m; 
        }
};
void CubicPPSplinescomplete::CalculateCoefficients(){
    int n=x.size();
    //使用eigen求解Am=b,（lamda_i*m_i+2m_i+mu_i*m_{i+1}=3mu_i*f[x_i,x_{i+1}]+3lamda_i*f[x_{i-1},x_{i}]）
    //其中，lamda_i=(x[i+1]-x[i])/(x[i+1]-x[i-1]),mu_i=(x[i]-x[i-1])/(x[i+1]-x[i-1]),和为1；
    std::vector<double> lamda(n-2);
    std::vector<double> mu(n-2);
    for(int i=0;i<n-2;i++){
        lamda[i]=(x[i+2]-x[i+1])/(x[i+2]-x[i]);
        mu[i]=(x[i+1]-x[i])/(x[i+2]-x[i]);
    }
    Eigen::MatrixXd A=Eigen::MatrixXd::Zero(n,n);
    A(0,0)=1;
    A(n-1,n-1)=1;
    for (int i=1;i<n-1;i++){
        for(int j=1;j<n-1;j++){
            if(j==i-1){
                A(i,j)=lamda[i-1];
            }else if(j==i){
                A(i,j)=2;
            }else if(j==i+1){
                A(i,j)=mu[i-1];
            }
        }
    }
    Eigen::VectorXd b=Eigen::VectorXd::Zero(n);
    b(0) = m_1;
    b(n-1) = m_n;
    b(1)=3*mu[0]*(f[2]-f[1])/(x[2]-x[1])+3*lamda[0]*(f[2]-f[1])/(x[2]-x[1])-lamda[0]*m_1;
    b(n-2)=3*mu[n-3]*(f[n-1]-f[n-2])/(x[n-1]-x[n-2])+3*lamda[n-3]*(f[n-1]-f[n-2])/(x[n]-x[n-1])-mu[n-3]*m_n;
    for (int i=2;i<n-2;i++){
        double b_1=3*mu[i-1]*(f[i+1]-f[i])/(x[i+1]-x[i]);
        double b_2=3*lamda[i-1]*(f[i]-f[i-1])/(x[i]-x[i-1]);
        b(i)=b_1+b_2;
    }
    //计算m矩阵
    Eigen::VectorXd m_A=Eigen::VectorXd::Zero(n);
    m_A = A.colPivHouseholderQr().solve(b);
    //存储m
    m.resize(n);
    for (int i = 0; i < n; i++)
    {
        m[i]=m_A(i);

    }
    //s'(a)=f'(a),s'(b)=f'(b)
    poly.resize(n-1);
for(int i=0; i<n-1; i++){
    double K_i=(f[i+1]-f[i])/(x[i+1]-x[i]);   
    Polynomial polys;
    polys.left_point = x[i];
    polys.right_point = x[i+1];
    polys.degree = 3;
    // 调整coefficients大小
    polys.coefficients.resize(4);  // 因为是3次多项式，需要4个系数
    
    // 或者使用push_back
    polys.coefficients.clear();  // 清除之前的数据
    polys.coefficients.push_back(f[i]);  
    polys.coefficients.push_back(m[i]);
    polys.coefficients.push_back((3*K_i-2*m[i]-m[i+1])/(x[i+1]-x[i])); 
    polys.coefficients.push_back((m[i]+m[i+1]-2*K_i)/pow(x[i+1]-x[i],2));
    poly[i] = polys; 
}

}


//自然PP样条
class CubicPPSplinesnatural:public PPSplines{
    private:
    double m_1;
    double m_n;
    double M_1=0;//左端点二阶导
    double M_n=0;//右端点二阶导
    std::vector<double> M;
    public:
    CubicPPSplinesnatural(const std::vector<double>& x,const std::vector<double>& f,const double& m_1,const double& m_n):PPSplines(x,f),m_1(m_1),m_n(m_n){
        CalculateCoefficients();
    }
    void CalculateCoefficients() override;
};

void CubicPPSplinesnatural::CalculateCoefficients(){
    int n=x.size();
    std::vector<double> lamda(n-2);
    std::vector<double> mu(n-2);
    for(int i=0;i<n-2;i++){
        lamda[i]=(x[i+2]-x[i+1])/(x[i+2]-x[i]);
        mu[i]=(x[i+1]-x[i])/(x[i+2]-x[i]);
    }
    //使用eigen求解Am=b,（lamda_i*M_{i+1}+2M_i+mu_i*M_{i-1}=6f[x_{i-1},x_i,x_{i+1}]）
    //其中，lamda_i=(x[i+1]-x[i])/(x[i+1]-x[i-1]),mu_i=(x[i]-x[i-1])/(x[i+1]-x[i-1]),和为1；
    Eigen::MatrixXd A=Eigen::MatrixXd::Zero(n,n);
    A(0,0)=1;
    A(n-1,n-1)=1;
    for (int i=1;i<n-1;i++){
        for(int j=1;j<n-1;j++){
            if(j==i-1){
                A(i,j)=mu[i-1];
            }else if(j==i){
                A(i,j)=2;
            }else if(j==i+1){
                A(i,j)=lamda[i-1];
            }
        }
    } 
    Eigen::VectorXd b=Eigen::VectorXd::Zero(n);
    b(0)=0;
    b(n-1)=0;
    for(int i=1;i<n-1;i++){
        double m_i=(f[i]-f[i-1])/(x[i]-x[i-1]);
        double m_i1=(f[i+1]-f[i])/(x[i+1]-x[i]);
        b[i]=6*(m_i1-m_i)/(x[i+1]-x[i-1]);
    }
    //计算M矩阵
    Eigen::VectorXd M_A=Eigen::VectorXd::Zero(n);
    M_A=A.colPivHouseholderQr().solve(b);
    //存储M
    M.resize(n);
    for(int i=0;i<n;i++)
    {
        M[i]=M_A(i);
    }
    CubicPPSplinescomplete complete(x,f,m_1,m_n);
    std::vector<double> m=complete.getFirstDerivatives();
    poly.resize(n-1);
    for(int i = 0; i < n-1; i++) {
    Polynomial polys;
    polys.left_point = x[i];
    polys.right_point = x[i+1];
    polys.degree = 3;
    polys.coefficients.clear();  // 清除之前的数据
    polys.coefficients.push_back(f[i]);
    polys.coefficients.push_back(m[i]);
    polys.coefficients.push_back(M[i]/2);
    polys.coefficients.push_back((M[i+1]-M[i])/(6.0*(x[i+1]-x[i])));
    poly[i] = polys;  // 使用索引赋值
}
}
//周期PP样条
class CubicPPSplinesperiodic : public PPSplines {
private:
    std::vector<double> M; // 二阶导数
    std::vector<double> m; // 一阶导数

public:
    CubicPPSplinesperiodic(const std::vector<double>& x, const std::vector<double>& f) : PPSplines(x, f) {
        CalculateCoefficients();
    }
    void CalculateCoefficients() override;
};

void CubicPPSplinesperiodic::CalculateCoefficients() {
    int n = x.size();

    // 检查周期性条件
    if (std::abs(f[0] - f[n - 1]) > 1e-8) {
        std::cout << "f(a) must equal f(b)!" << std::endl;
        return;
    }

    // 计算 lamda 和 mu
    std::vector<double> lamda(n - 2);
    std::vector<double> mu(n - 2);
    for (int i = 0; i < n - 2; i++) {
        lamda[i] = (x[i + 2] - x[i + 1]) / (x[i + 2] - x[i]);
        mu[i] = (x[i + 1] - x[i]) / (x[i + 2] - x[i]);
    }

    // 构造矩阵 A 和向量 b
    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(n, n);
    Eigen::VectorXd b = Eigen::VectorXd::Zero(n);

    // 周期性边界条件
    A(0, 0) = 1;
    A(0, n - 1) = -1; // f'(a) = f'(b)
    A(n - 1, 0) = 1;
    A(n - 1, n - 2) = -1; // f''(a) = f''(b)

    // 填充 A 的中间部分
    for (int i = 1; i < n - 1; i++) {
        A(i, i - 1) = mu[i - 1];
        A(i, i) = 2;
        A(i, i + 1) = lamda[i - 1];
    }

    // 填充 b 的中间部分
    for (int i = 1; i < n - 1; i++) {
        double m_i = (f[i] - f[i - 1]) / (x[i] - x[i - 1]);
        double m_i1 = (f[i + 1] - f[i]) / (x[i + 1] - x[i]);
        b(i) = 6 * (m_i1 - m_i) / (x[i + 1] - x[i - 1]);
    }

    // 计算二阶导数 M
    Eigen::VectorXd M_A = A.colPivHouseholderQr().solve(b);
    M.resize(n);
    for (int i = 0; i < n; i++) {
        M[i] = M_A(i);
    }

    // 计算一阶导数 m
    m.resize(n);
    for (int i = 0; i < n - 1; i++) {
        m[i] = (f[i + 1] - f[i]) / (x[i + 1] - x[i]) - (M[i + 1] + 2 * M[i]) * (x[i + 1] - x[i]) / 6.0;
    }
    m[n - 1] = m[0]; // 周期性条件

    poly.resize(n-1);
    for(int i = 0; i < n-1; i++) {
    Polynomial polys;
    polys.left_point = x[i];
    polys.right_point = x[i+1];
    polys.degree = 3;
    polys.coefficients.clear();  // 清除之前的数据
    polys.coefficients.push_back(f[i]);
    polys.coefficients.push_back(m[i]);
    polys.coefficients.push_back(M[i]/2);
    polys.coefficients.push_back((M[i+1]-M[i])/(6.0*(x[i+1]-x[i])));
    poly[i] = polys;  // 使用索引赋值
}
}
#endif
